/* Geometry.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#if !defined _geometry_
#define _geometry_
/*! An essay to define geometrical object and distance to them.
*/

/*! A factory mechanism to crete geometric object.
*/
template<class data>
struct geometric_object : public boost::bitwise<geometric_object<data> >
{
	typedef geometric_object type;
	typedef geometric_object Base_Cast;
	typedef smart_ptr<Base_Cast> smart_ptr_type;
	type(){};
	type(const type &from) {
			first_ = from.Create_Instance();
	};
	type& operator = (const type &from) {
		if (from.first_)
				first_ = from.first_->Create_Instance();
		return *this;
	};
	template <class T>
	smart_ptr_type Create_Instance_helper() const
	{
		T* result = new T;
		if (first_)
			result->first_ = first_->Create_Instance();
		return smart_ptr_type(result);
	};
	template <class T, class Data>
	smart_ptr_type Create_Instance_helper_arg(const Data * data) const
	{
		smart_ptr_type result(new T(*data));
		if (first_)
			result->first_ = first_->Create_Instance();
		return result;
	};



	virtual double distance(const data & point) const {
		return first_->distance(point);
	};
	virtual data grad_distance(const data & point) const {
		return first_->grad_distance(point);
	};
	virtual smart_ptr_type Create_Instance() const {
		if (first_)
			return first_->Create_Instance();
		return smart_ptr_type();
	};
	type& operator |=(const type & from) {
		if (first_)
			first_ = smart_ptr_type(new OS_or<data>(*first_,from));
		else first_ = from.Create_Instance();
		return *this;
	};
	smart_ptr_type first_;
};

/*! A tool to compute binary operation over geometric object
*/
template<class data>
struct OS_Geometry_binary_operator : public geometric_object<data>
	{
		typedef geometric_object				Base;
		typedef Base::Base_Cast					Base_Cast;
		OS_Geometry_binary_operator(){};
		OS_Geometry_binary_operator(const Base &first) : Base(first)
		{};
		OS_Geometry_binary_operator(const Base &first, const Base &second) {
			first_ = first.Create_Instance();
			second_ = second.Create_Instance();
		};

		template <class T>
		smart_ptr_type Create_Instance_helper() const
		{
				T* result = new T;
				if (second_)
					result->second_ = second_->Create_Instance();
				if (first_)
					result->first_ = first_->Create_Instance();
				return smart_ptr_type(result);
		};
		smart_ptr_type second_;		
	};
/*! Union of geometrical objects
*/
template<class data>
	struct OS_or
		: public OS_Geometry_binary_operator<data>
	{	
		typedef OS_Geometry_binary_operator	Base;
		OS_or(){};
		OS_or(const Base_Cast &first, const Base_Cast &second) 
			: Base(first,second)
		{};
		virtual std::string get_name() const {return "|";};
		virtual double distance(const data & point) const {
			return std::min( first_->distance(point),  second_->distance(point) );
		};
		virtual data grad_distance(const data & point) const {
			double d_1 = first_->distance(point);
			double d_2 = second_->distance(point);
			if (d_1<d_2)
				return first_->grad_distance(point);
			return second_->grad_distance(point);
		};

		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper<OS_or>();
		};
	};

/*! A circle object
*/
template <class data>
struct Circle : public geometric_object<data>
{
	typedef Circle type;
	typedef geometric_object Base;
	Circle(){};
	Circle(const Circle & from) : center_(from.center_), radius_(from.radius_){};
	Circle(const data center, const OS_double radius) : center_(center), radius_(radius)
	{};
	OS_double distance(const data& from) const
	{
		double norm = jmmath::norm2(from-center_);
		return std::abs(norm-radius_);
	};
	data grad_distance(const data& from) const
	{
		data result(from-center_);
		double norm = jmmath::norm2(result);
		if (norm == 0.) {
			result.assign(std::sqrt(1./result.size() ) );
			return result;
		}
		result /= norm;
		if (norm == radius_)
			return result;
		result *= ( (norm -radius_) / std::abs(norm -radius_));	
		return result;
	};
	virtual smart_ptr_type Create_Instance() const {
		return Create_Instance_helper_arg<type,type>( this ); 
	};
	data center_;
	OS_double radius_;
};

/*! A segment object
*/
template <class data>
struct Segment : public geometric_object<data>
{
	typedef Segment type;
	typedef geometric_object Base;
	Segment(){};
	Segment(const Segment & from) : left_(from.left_), right_(from.right_), left_right_(from.left_right_), left_right_norm_(from.left_right_norm_){};
	Segment(const data& left, const const data& right) : left_(left), right_(right)
	{
		init();
	};
	void init() {
		left_right_ = left_-right_;
		left_right_norm_ = jmmath::norm2(left_right_);
		left_right_ /= jmmath::norm2(left_right_);
	};
	OS_double distance(const data& from) const
	{
		double scalar = (from-right_)*(left_right_);
		if (scalar <= 0)
			return jmmath::norm2(from-right_);
		if (scalar >= left_right_norm_)
			return jmmath::norm2(from-left_);

		return jmmath::norm2((from-right_)-scalar*left_right_); 
	};
	data grad_distance(const data& from) const
	{
		double scalar = (from-right_)*(left_right_);
		if (scalar < 0)
			return (from-right_)/(jmmath::norm2(from-right_));
		if (scalar > left_right_norm_)
			return (from-left_)/(jmmath::norm2(from-left_));

		data B = (from-right_)-scalar*left_right_;
		scalar = jmmath::norm2(B);
		if (scalar != 0.)
			return B / jmmath::norm2(B);
		B.assign(0.);
		return B;
	};
	virtual smart_ptr_type Create_Instance() const {
		return Create_Instance_helper_arg<type,type>( this ); 
	};
	data left_;
	data right_;
	data left_right_;
	double left_right_norm_;
};


#endif